function [cInt,cArea] = extractMetrics(I,BW,channels,positions,time)

%%Constants
area_min = 20; %single cell
area_max = 512*512/3; %flood

%% Intensities
if nargin < 5
    time = 1:size(I{1,1},3);
    if nargin < 4
        positions = 1:size(I,2);
        if nargin <3
           channels = [2 3]; 
        end
    end
end

h_bar = waitbar(0,'Quantifying Features...');

for p=positions
    for t=time
        %Cluster masks from combined fluorescent channels
        cc = bwconncomp(BW{2,p}(:,:,t)|BW{3,p}(:,:,t),8); %8-connected
        lcc = labelmatrix(cc);
        cellstats = regionprops(cc,'Area');
        area = cell2mat({cellstats.Area}');
        %Set min area cutoff
        area(area<area_min)=[];
        %Set max area cutoff (e.g. flooded frame)
        area(area>area_max)=[];
        for i=1:numel(area)
            for c=2:3
                cPx{c,p,i}(:,:,t) = I{c,p}(:,:,t).*uint16(lcc==i);
                cInt(c,p,i,t) = sum(sum(cPx{c,p,i}(:,:,t)));
                cArea(c,p,i,t) = area(i); 
            end
        end
        %Status
        waitbar(sum([p-1 t].*[numel(time) 1])/prod([numel(positions) numel(time)]),h_bar)
        clear cc lcc
    end
end
cInt(cInt==0)=NaN;
cArea(cArea==0)=NaN;
close(h_bar)

end